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Abstract 

We suggest to include the density of electron charge explicitly in the 
electron potential of density functional theory, rather than implicitly via 
exchange-correlation functionals. The advantages of the approach are 
conceptual and numerical. Conceptually, it allows to formulate a physical 
principle for the development of quantum mechanical systems: it is the 
principle of energetic equilibrium, because the energy principle, in this 
case, applies not only globally, but also on a local level. The method is an 
order- TV method, scaling linearly with the number of atoms. It is used to 
calculate the electronic groundstate of a metallic surface, where we find 
good agreement with experimental values. 



1 Density functional theory 

A key innovation in theoretical solid state physics in the last fifty years was 
the reformulation of quantum mechanics in a density formalism, based on the 
Hohcnbcrg-Kohn theorem pp. Despite initial resistance, in particular from quan- 
tum chemists, the method has replaced previous frameworks and provides, to 
date, the most advanced theoretical model for the calculation of atoms, solids, 
and liquids. However, its implementation relies on a rather cumbersome de- 
tour. While the Hohenberg-Kohn theorem is formulated exclusively in terms of 
electron charge densities and energy functionals, calculations today are based 
almost exclusively on the specifications given by Kohn and Sham one year after 
the initial theorem was made public PJ. These Kohn-Sham equations define, 
how the electron interactions within a system can be split into three separate 
single-electron parts: a Coulomb potential Ucom giving the electron-ion attrac- 
tion and the electron-electron repulsion; an exchange potential Ux, taking into 
account the Pauli exclusion principle; and a correlation potential Uc, covering 



the correlated motion of electrons in a many-electron system. Then every elec- 
tron in the system is described by a single-particle Schrodinger equation with 
the effective potential U e ff, 



U ef f =U C ou + Ux +U C 



(1) 



in order to determine the charge density. Within the local density approxima- 
tion, the exchange and correlation potentials are usually combined and thought 
to be depending on the local density of electron charge p(r): 



The variational problem thus contained at every single step of the energy mini- 
mization cycle a solution of the single-particle Schrodinger equation, from which 
the electron eigenstates and their eigenvalues are determined. The density of 
charge at a given step of the cycle is calculated by adding single-electron charges: 



N denotes the total number of electrons in the system. While this procedure 
is generally successful, and implemented today in numerical methods optimized 
for efficiency (see e.g the ingenious way ionic and electron degrees of freedom 
are treated on much the same footing following a method developed by Car and 
Parinello jH]), it is highly inefficient in one crucial conceptual point: If, according 
to the Hohenberg-Kohn theorem, the density of charge is the only physically 
relevant variable of the system, then solving the Schrodinger equation, setting up 
the eigenvectors, and computing the density of electron charge is an operation, 
which creates a vast amount of redundant information. Every information, 
pertaining to the solution of the single-particle Schrodinger equation and the 
summation of single electron charges is discarded at the end of every step in the 
iteration cycle. One could therefore say that more than 90% of the information 
created in today's simulations is actually irrelevant. The question thus arises: 
Do we have to create this information at all, or can we find a more direct way 
to arrive at the groundstate density of charge without this cumbersome detour 
via the single-particle Schrodinger equation? 

From a physical point of view, this state of affairs is easy to understand, if one 
considers the basic innovation, quantum mechanics introduced into the original 
framework of mechanics. It is the property of phases and phase correlations, 
which, even though it is present in any many-electron system, does not find an 
expression in its key quantity, the density of charge. Inducing phase correlations 
into density functional theory thus meant to rebuild a phase coherent system 
by way of the original Schrodinger equation. The easiest way to build such a 
system is to determine the coherent wavefunctions of single electrons and to 
ensure that the states are orthogonal in Hilbert space. 

A more general approach would be to look at phase-correlations within a 
different framework. This framework, the theory of density matrices @JE1E]> 



Uxc(r)=U xc (p(r)) 



(2) 



N 




(3) 
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has consequently been the focus of research over a period of more than ten 
years. However, a framework combining the notion of phase-coherence for 
many-electron wavef unctions, and the basic requirement in DFT to find the 
groundstate density, has so far been missing. In this paper, we introduce such 
a direct method and show how it can be used to determine the groundstate 
density without any recourse to the single-particle Schrodinger equation. 

The paper was written in the spirit of orbital free density functional theory 
003 EH- It is frequently pointed out by protagonists of this line of research, that 
an orbital free, and Bloch-state free formulation of DFT has enormous potentials 
for numerical improvements. The main difference to the present framework is 
philosophical, rather than practical. While orbital free developments are essen- 
tially based on the Hohenberg-Kohn theorem and attempt to find a transferable 
formulation for the kinetic energy functional, we have analyzed the foundation 
of the theory itself and sought a continuous reformulation of the fundamen- 
tal theory. Not least, because such a reformulation has substantial conceptual 
advantages, as will be shown in the following sections. 

The paper is organised as follows: in section II we introduce the theoretical 
concept, based on a continuous generalisation of potentials in the many-electron 
Schrodinger equation. It is shown that the ensuing equation for the density of 
charge is non-linear, and that Schodinger's equation is a zero-order approxima- 
tion of the general relation. In section III we sketch the solution of the eigenvalue 
problem for a solid-state system. Finally, the new formulation is tested for metal 
surfaces in section IV. 

2 Theoretical basis 

The theorem of Hohenberg and Kohn elevated the density of charge to the state 
of a physical variable in many-electron systems. In addition, observations by 
scanning tunnelling microscopes seem to confirm that the observed structures 
at a surface are essentially due to the distribution of electron charge, measured 
with a resolution of a few picometer ^UJ- It is then quite astonishing to find that 
the electron charge has to be determined by a complicated computation process 
from the single electron states. We show in the following, how a representation 
of the system, without any reference to these states, can be constructed on the 
basis of Green's functions. 

2.1 Green's functions and phase correlations 

The Green's function of a many-electron system can be expanded into single 
electron states by a spectral decomposition ^] ^] : 



Here, G(r, r',z) is the many-electron Green's function of the system in a 
local basis, N the number of electrons, and z the complex energy parameter. 
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The Green's function will become its complex conjugate upon an exchange of 
the local coordinates: 

G{v,v',z) = G*{v',v,z*) (5) 
It complies with a differential equation of the Hamiltonian in a local basis: 

(H(r, r') - z) G(r, r', z) = -S(r - r') (6) 

Note that the Green's function of the many-electron system depends only on 
two local coordinates r and r'. We show now that this Green's function can be 
written as the product of the phase correlation functions \l/(r) and ^(r') and 
an energy function F(z): 

G(r,r',z) = F(z)y(r)V*(r') (7) 

where |\I/(r)| 2 is equal to the density of charge upon a suitable choice of F(z). 
First, it follows from Eq. [7| that the conjugate complex of the F(z)5 , (r)\E'*(r') 
complies with the same relation as the Green's function, i.e. it is equal to the 
transposed Green's function: 

G*(r', r, z*) = [F(2*)#(r')#*(r)]* = F(z)**(r')*(r) = G(r, r',z) (8) 

Taking the diagonal Green's function in the limit of real eigenvalues gives: 

G ± (t,t,E) = *(r)**(r) lim o F(E ± irf) (9) 

The imaginary part of the diagonal Green's function, integrated over energy, 
must be equal to the charge density: 

1 

p(r) = T- / dEImG ± {v,v,E) 

1 r + °° 

--- / dEIm lim F{E ± ir?)*(r)**(r) (10) 
7T J_oo ? '^+° 



=F- 



If ^'(r)^'*(r) = p(r), the rest on the right side must be unity. This is satisfied 
if 

T-Im lim F{E±irf) = S(E - Eq) (11) 

7T t?^+0 



Using the fact that 



y^+o x ±iy 
gives us the explicit form of F(z 



lim — — = P ( - ) =F iir5{x) (12) 



m = —pr (is) 

Z — t,n 
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Thus, the Green's function can be written as the product of two local phase 
correlation functions [12]: 

Z — £/q 

At this stage we have shown that the phase correlation function ^(r) is an 
exact representation of the TV-electron system, as long as the system can be 
described by a Green's function Eq. 0] However, it is not yet clear, whether 
there is a differential equation describing the evolution of ^(r). This equation 
cannot be equal to the standard Kohn-Sham equations, as this would lead to a 
non-local concept. If we write the Hamiltonian H(r,r') as the sum: 

H(r,r') = H(r) + H(r') - ff 1 (r)irj(r / ), (15) 

and insert into Eq. [5] we get: 

**(r') [H(r) - z] *(r) *(r) [H(r') - z] **(r') 



z - E z - E 

[Hi(r)HZ(r') - z]*(r)**(r') 



(16) 



En 



-S(r-r' 



It can be seen that the second line contains a product of operators at different 
locations r and r'. This would mean that the phase correlation at a point 
r depends on events at point r', even though there is no obvious connection. 
However, since ^ depends only on r, this seems unjustified. In addition, the 
Hohenberg-Kohn theorem proves that if ^(r)! 2 is equal to the groundstate 
density of charge, then an energy eigenvalue Eq is defined and must be equal 
to the total energy. Therefore it should be possible to formulate the problem of 
finding the phase correlation function 'J'(r) in an eigenvalue equation. In this 
case the many-electron Green's function has only one pole, in contrast to the 
standard spectral decomposition, where it has N poles. 

Here we make the only conjecture in this paper: we assume that the phase 
correlation function \&(r) is described by a suitable modification of a single 
coordinate Schrodinger equation: 

_|_V 2 *(r) + C/(r)*(r) = E V(r) (17) 
2m 

with a local potential U(r). That such a modification is possible and leads 
to correct results in simple test cases, will be shown in the following. The 
conjecture is supported by the derivation of an approximate single coordinate 
Schrodinger equation from the Hartree-Fock model of an TV-electron system |13j . 

If this conjecture is justified, then we have obtained, at this point, a phase 
correlation function, which includes all phase correlations in the occupied range, 
it complies with a local equation, its corresponding eigenvalue Eq is equal to the 
total energy of the system, and its square gives the electron charge density in 
the system. This entity looks very much like a descriptor of the physical system 
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itself. Something, which is otherwise called a many-body wavefunction. As this 
term might lead to confusion, since the many-body wavefunction is in standard 
theory described by N local coordinates (ri, r 2 , r 3 , rjv), we shall call it the 
many-electron wavefunction in the rest of the paper. 

An additional note concerning the relation of the framework presented in 
this paper to density functional theory seems necessary. In particular, since 
the original derivation of the Hohenberg-Kohn theorem is based on many-body 
wavefunctions [Q. In this respect we note that the groundstate charge density, 
once obtained, can be written in any representation. Thus an energy minimum, 
associated with a density, is enough to describe the physical content of a many- 
electron system. As long as the energy is a minimum, it is therefore irrelevant, 
how this minimum was obtained in practice. This is reflected by the most 
efficient methods in density functional theory, which are based on trial wave- 
functions, which do not diagonalize the Hamiltonian. Instead, the off-diagonal 
elements are minimized with the help of numerical algorithms until they vanish. 



2.2 The non-linear Schrodinger equation 

So far, we have not considered the modified form of the Hamiltonian H(r). In 
particular the potential, which in the standard formulation is the potential of 
single electrons, is not suitable for our purposes. To find a local formulation 
of the problem, we start with the Schrodinger equation of the many-electron 
system, written as: 

H 2 

V 2 *(r) + [/(r)*(r) = EV(r) (18) 

Here, V 2 denotes the Laplacian, U(r) a potential in the system, and E the total 
energy. Within the context of density functional theory it is generally assumed 
that the potential U reflects the interaction between the charge of one electron, 
and the potential of the environment. This entails that the potential is a sum 
containing the positions of all electrons in the system. While such an approach 
is possible, it does not reflect the previous findings. If phase correlations and 
charge are local quantities, then the interactions between electron charge and 
the potential of all other charges should also be local. 

The single coordinate equation 1181 must therefore be recast into a relation 
based on densities, and not integral quantities. This provides the desired equa- 
tion for the many-electron wavefunction. If U(r) is e.g. the Coulomb potential 
of an ion with charge Z , then the potential density of electron attraction at a 
point (r) is given by: 

, x , 1 Ze 2 p(r) , , 1 Ze 2 
v(r)p(r) = ■. v(r) = ■. r (19) 

Similarly, the kinetic energy term and the total energy will change into kinetic 
energy density and total energy density. The variables in the local equation will 



G 



be changed into: 



— =>2f U(r) v(t) P (t) E^e (20) 

/ is a system parameter, related to the volume of the system, which shall be 
analysed in detail further down. v(r) is the external potential, e is the energy 
density. The equation then reads: 

-2/V 2 *(r) + w(r)p(r)tf (r) = etf (r) (21) 



2.2.1 External potential 

The interaction between electrons and ions in a groundstate many-body problem 
is exclusively described by Coulomb interactions. If \P(r), as shown in previous 
sections, uniquely describes the problem, then the same should hold for the 
modified Schrodinger equation. Exchange correlation potentials, as in density 
functional theory, would enter the description only at the level of single electron 
states. Since ^(r) describes the total system, and not single electron state, they 
will not enter the present framework. However, Coulomb interactions will be 
modified, if one considers magnetic systems. In this case the potential v(r) will 
depend on the spin-state of electron charge. At present, we do not consider 
this case. The potential v(r) then is the Coulomb interaction potential due to 
electrons and ions, 



w(r) 



Atteq 



|r-r'| £f|r-Ri 



(22) 



2.2.2 Charge density 

The density of charge p(r) is the square of the many-electron wavefunction, or: 



P (t) = **(r)*(r) (23) 

Writing an equivalent equation for the complex conjugate function \P*(r), we 
get: 

-2/V 2 tf*(r) + v(r)p(r)tf*(r) = e**(r) (24) 

Multiplying the first of these equations by \£*(r), and the second by ^(r), and 
adding the equations, we end up with: 

-/ [tf*(r)V 2 tf (r) + *(r)V 2 **(r)] + v(r)p(r)p(r) = ep(r) (25) 

The Laplacian acting on the density is given by: 

V 2 [*(r)**(r)] - [**(r)V 2 *(r) + *(r)V 2 **(r)] + 
+ 2[V*(r)] [V**(r)] (26) 

The second term on the right contains the essential phase correlations of the 
many-electron wavefunction. We symbolize it by II[*] to denote that it contains 



7 



the many-electron phase correlations. In principle there seems no straightfor- 
ward method to describe it in terms of charge densities alone: as the many- 
electron system has to be phase coherent, the density of charge alone provides 
insufficient information. If this were not the case, then a formulation based on 
density alone would be sufficient to describe the physical content of a many- 
electron system. From this viewpoint, the essential finding of density functional 
theory is not only the uniqueness of the charge distribution, but also of the 
groundstate energy. This energy, in turn, can only be determined by construct- 
ing a phase coherent system. Thus the density of charge will be described by: 

[-/V 2 + W (r)p(r)]p(r) + /n[vI/] = ep(r) 

= 2[V*(r)][V**(r)] (27) 

However, the non-linear equation describes the many-electron wavefunction, 
provided the potential and groundstate charge density are known. It is therefore 
not strictly necessary to find a formulation based on density alone, as long as the 
nonlinear equation leads to correct many-electron wavefunctions. The derivation 
thus provides a set of three equations, defining the groundstate of a system: 

[-2/V 2 +v(r)p(r)] *(r) = e#(r) 
[-/V 2 + Kr)p(r)]p(r) + /n[*] = ep(r) (28) 

p(r) = *(r)**(r) 

Note that this set of equations applies to the many-electron density of charge 
without any reference to single Kohn-Sham states. It describes the physical sys- 
tem by four variables: the density of charge p(r), the external potential v(r), 
which will also depend on the density, the energy density e, and the many- 
electron wavefunction \&(r). The equations lead to a direct method of cal- 
culating the groundstate charge density by solving the nonlinear Schrodingcr 
equation and minimizing the energy density e of the system. The self-consistent 
procedure, elaborated in the following section, will be 

• make an initial guess for the density p , 

• construct the potential vq, 

• solve the nonlinear Schrodingcr equation for and 

• construct a new charge density p\ = \I r * 1 ^^ 

The cycle is repeated until input and output charge densities are equal. The 
problem of finding the groundstate energy density can be formulated as the 
following minimization problem: 

_ f J rf 3 r ( [-/V 2 + v(r)p(r)} p(r) + /n[tt(r)]) 

Note that also in this case the problem cannot be formulated as a minimization 
problem of charge density alone. This is a consequence of basing the whole 
derivation on the many-electron problem and ensuring the phase coherence 
throughout the system. 
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2.3 Physical implications 

From a dimensional point of view it is interesting to analyse the meaning of the 
system constant /. Its dimension is that of a force: 



[/v 2 ] = 



E 

1? 



[/] = 



E 



N(SI) (30) 



Apart from the universal constants h and m the system constant / only depends 
on the volume. In this respect, it is beyond the well known Thomas-Fermi model, 
where the kinetic energy density is system independent. The second difference 
is that phase correlations arc included in the picture, by solving the non-linear 
Schrodinger equation. The differences are crucial. As shown in the following 
sections, the model recaptures the density oscillations at metal surfaces, con- 
trary to the Thomas-Fermi model. That the kinetic energy density cannot be 
system independent, is in line with observations in quantum mechanics that the 
energy of a system comprising e.g. a single electron is inverse proportional to 
the enclosing volume. It seems therefore that the action of an external potential 
on the density of charge can be seen as a force acting on the wavevector, or the 
spatial distribution of the density. The reason for this force is that energy at 
a given point of the system has to be conserved. The reason that it acts on 
the curvature rather than the amplitude is the relation between curvature or 
wavevector and kinetic energy, predicted by Louis de Broglie, and experimen- 
tally confirmed by Davisson and Germer |14l I15| . 

Since e was defined as the total energy divided by the volume the equation 
has an additional physical content. A solution of the equation with constant 
e means that the energy density does not vary from one point of the system 
to another. From the viewpoint of thermodynamics, this is also the condition 
of a system in equilibrium. One could thus say that the specific form of the 
Schrodinger equation has one very specific physical meaning: it provides the so- 
lution of distributing electron charge in an external potential in such a way, that 
the energy density is constant. It explains, why the mathematical form of the 
equation is different from e.g. a differential equation in classical electrodynam- 
ics. There, the boundary conditions usually lead to a selection of amplitudes 
and frequencies in a system. Here, the frequencies need to be modified at every 
single point in order to achieve constant energy density values. In principle, this 
is still a wave-dynamical problem, but a problem modified by the ability of elec- 
trons to interact with external potentials and with each other. This formulation 
of the many-electron problem also explains, why wave mechanics, based on the 
Schrodinger equation, appears to be non-local |16| : since the main condition 
implicitly contained in the equation is the condition of thermal equilibrium, one 
cannot separate one part of the system from any other part, without contra- 
dicting the assumption that the system as a whole will equilibrate. Thus, wave 
mechanics is in fact a local theory, like every other theory known in physics, 
but it is applied to systems which are assumed to be connected throughout. 
This can clearly be seen in the mathematical form of Eqs. |2H1 a h variables of 
the system, v(r),p(r), and e are local scalars, defined at a specific point of the 
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system. These local variables cannot, however, be arbitrarily chosen, since e 
must be a constant and / must relate to the volume of the system. 

One might object, at this point, that the problem is now vastly more com- 
plicated. After all, the eigenvalue problem is now non-linear, and it is generally 
impossible to compute solutions of non-linear differential equations in full gen- 
erality. Here, we remember the main result of the Hohenberg-Kohn theorem: 
the groundstate charge density is the charge density, which minimizes the total 
energy and thus the energy density e. This shall provide, as will be seen, a 
way to compute the exact solution starting from an initial guess for the charge 
density p(r), by solving a linear eigenvalue equation and adjusting the charge 
density from one iteration to the next. In principle, it is a very similar method 
to the one used today in standard DFT codes, but for two decisive advantages: 

1. The charge density is computed directly without any reference to single- 
particle Schrodinger equations. 

2. The charge density computed is the many-electron charge density and it is 
therefore completely independent of exchange-correlation functionals and 
density approximations. 

2.4 Possible objections 

From a pragmatic point of view, the approach seems not too different to the 
standard approach, apart from the fact that the non-local density-correlations 
in the potential have been replaced by the charge density itself. However, there 
is a number of issues, related to non-linear concepts, which deserve a closer look. 

• If the Hohcnberg-Kohn theorem is correct, then the groundstate charge 
density is unique. Here, we have a non-linear potential, and the ground- 
state is therefore potentially not unique. 

There are two possible answers. The first is the current state of DFT itself. After 
all, exchange correlation functionals include the density and its derivatives. The 
effective potential in today's DFT simulations is therefore also non-linear. This, 
it seems, does not invalidate the results of calculations. The second is based 
on recent research (in Mathematics) on the behavior of non-linear Schrodinger 
equations, related to weak coupling in a many-electron systems ^3]: m this 
case the equations still yield physically acceptable solutions, even though the 
existence of a solution depends on reasonable physical conditions. 

• The Thomas- Fermi model also attempts to calculate charge densities from 
a non-linear potential and a kinetic energy density. The same seems to 
be tried here. But we know that this model fails even to explain chemical 
bonding. 

The answer to this question lies in the treatment of the kinetic energy term, 
and the use of the nonlinear Schrodinger equation to determine many-electron 
charge densities. Within the Thomas-Fermi model, the prefactor of the kinetic 
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energy functional is a constant. This means, that it does not depend on the 
system, but is universally the same. Physically speaking, this cannot be correct. 
Because an electron in a very large negative potential must experience different 
kinetic effects than an electron in a comparatively small negative potential. In 
the present model the kinetic energy density is inversely proportional to the 
volume, and contains a term which is directly related to the many-electron 
phase correlations. It will be seen in the simulation of metal surfaces, that this 
leads to the correct behavior of electrons at a boundary. 

• I do not see spin-statistics involved in this approach. In standard DFT 
this is part of the exchange-correlation functionals. Where is it in this 
formulation? 

Here, one has to note that spin-statistics are also not involved in the formulation 
of the non-magnetic many-electron Schrodingcr equation. They only arise at 
the level, where single electrons are thought of as individual entities. A many- 
electron wavefunction in a many-electron Schrodingcr equation is a unique scalar 
function of location and spin. In a non-magnetic formulation only of location. 
So far, we have not included spin in the framework. This must be done in field 
mediated manner and will be the topic of future work. 

• Non-linearity is nothing new. It is already present in the standard tech- 
niques, since self-consistency cycles are essentially non-linear. The non- 
linearity lies therefore in the variational principle, but not in the funda- 
mental equations. 

While we do not dispute that the variation and the ensuing self-consistency 
cycles are non-linear, the fundamental equations, the Kohn-Sham equations, 
are still linear. If this were not the case, then the groundstate of a system could 
not be described as a superposition of N single electron states. While this, we 
think, is correct for a mean-field theory, it is not suitable within a many-electron 
framework. In this sense the described approach extends the non-linearity to 
the fundamental equations themselves, with the effect, that superpositions are 
no longer possible. The groundstate calculated with this approach cannot be 
decomposed into single-electron states, or only, if one uses an approximation. 

Before sketching the solution of the eigenvalue problem, let us consider some 
limiting cases, where one can show that the linear Schrodingcr equation is in 
fact the zero order approximation of the general problem. 

2.5 Zero order approximations 

The main areas, where the Schrodinger equation is used today, are the theory 
of atoms, and solid state theory. In solid state theory the most frequently 
employed DFT codes use periodic boundary conditions and are generally applied 
to systems with translational symmetry. In quantum chemistry, DFT methods 
usually focus on atomic orbitals and their representation by suitable basis sets. 
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Let us consider initially an atomic system, where the density of charge can 
be closely mimicked by an exponential function. Thus, with radial symmetry: 

P( r ) = Po exp (-/er) (31) 

The expansion of the exponential function gives: 

/ nr (nr) 2 (rer) 3 \ . , 

P(r)=Po[l-j I + { -^- -^r + -j (32) 

The zero order approximation assumes that p(r) = p = constant. Inserting p 
into Eq. 1211 and multiplying by the volume V, we get: 

-2/VV 2 *(r) + v nl (r)p VV(r) = E9(t) (33) 

Setting 2fV = h 2 /2m we arrive at the conventional Schrodinger equation, e.g. 
for the hydrogen atom: 

*(r) = E#(r) (34) 

In solid state physics, systems are generally periodic. The density of charge in 
this case will be a periodic function. For the sake of demonstration we use a 
cosine, but the same argument applies to a sine or an imaginary exponent. Say, 
we can write the density of charge as: 

p(r) — a(l + /3cos (kr)) = a + 

where k is a vector in reciprocal space, and (3 < 1. We set a + a ■ j3 = pq. Then 
the zero order approximation for the non-linear eigenvalue problem will read 
(we have again multiplied by the volume V): 

-2/W 2 *(r) + v nl (r)p VV(r) = m(r), (36) 

and it will again lead to the conventional Schrodinger equation if we set 2fV = 
h 2 /2m: 

*(r) = E*(r) (37) 

We are well aware that these examples do not constitute a proof that the non- 
linear eigenvalue problem is exactly equivalent to the Schrodinger equation. In 
fact, if this were the case, then the whole formulation might be considered 
redundant, since it is only a different way of stating what is already known. 



_^V 2 + C/(r) 
2m 



_^V 2 + C/(r) 
Im 
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2.6 Removing the paradoxes from quantum mechanics 

The new framework removes two of the fundamental problems in quantum me- 
chanics, which have been the subject of debate almost from the beginning of 
modern physics. It is probably not exaggerated to say that the research done on 
quantum paradoxes has, over the years, produced a volume of work comparable 
to the volume of work in standard developments. 

1. The Schrodinger cat paradox ^7]. In the standard formulation the cat 
is either dead, or alive, but its state can only be determined by a mea- 
surement. There is, fundamentally, no objective way of deciding before a 
measurement. 

2. The Einstein-Podolsky-Rosen paradox JSj- I n this case a quantum me- 
chanical system is thought to be split into two subsystems which increase 
their distance with the velocity of light. After some time two separate 
measurements are performed on either subsystem, and according to the 
standard formulation, the measurement on one subsystem should deter- 
mine the measurement on the second one, even though no physical inter- 
action is known which exceeds the speed of light. 

Both paradoxes rest on one fundamental feature of systems in standard quan- 
tum mechanics: the first, on the possibility of superimposing separate solutions 
of the Schrodinger equation to describe the state of a system; the second, on the 
very same feature with a slight modification concerning the times of measure- 
ment and the critical distance between the subsystems. The resolution of the 
paradoxes, within this new framework, is the easiest possible answer: the cat 
is dead, or it is alive, independent from any measurement, because there exists 
at any given moment only one unique solution to the fundamental equations. 
This is a consequence of the non-linearity of the general framework, as already 
pointed out. The same applies to the second paradox, with the slight modifica- 
tion, that performing a measurement on one subsystem changes the state of this 
subsystem. While in the standard framework, this has an effect on the second 
subsystem, since the state of the overall system is described by a superposition, 
it does not affect the second subsystem in the present framework, since the fun- 
damental principle of equilibration would be violated. The second subsystem 
cannot equilibrate with the first one, because it cannot interact. 

From a physical point of view, this resolution of the paradoxes seems ac- 
tually more in keeping with a materialistic mindset. Physicists of a classical 
disposition have been arguing against the consequences of the mathematical 
framework, and in particular the feature of superposition, for almost a century. 
Here, we adopt the view that this consequence hinges on the linearity of the 
Schrodinger equation, and that this linearity is actually a feature enforced by 
the approximations used. 
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3 Iterative solution of the non-linear eigenvalue 
problem 

While some eighty years ago the non-linear problem would have been unsolv- 
able, this has changed significantly over the last fifty years due to the existence 
of powerful computers. In this respect, theory is not much different from exper- 
iments: a change in the equipment usually leads to a different focus in research, 
as phenomena, which were previously out of range, can now be analysed in a 
scientific manner. 

In particular the advent of DFT and the continuous improvement of com- 
puting algorithms applied to solid state physics have created an environment 
in theory, where the inability to solve equations analytically is no hindrance. 
Two of the main problems, needed to solve the non-linear problem, are already 
incorporated in standard DFT codes: (i) The problem of finding the charge 
density pi+\ after a given iteration i. This is done by elaborate mixing schemes, 
which are constructed to speed up the convergency process, (ii) The problem 
of finding the minimum energy density. This problem is not different from the 
problem of finding it in standard DFT convergency cycles, where one aims at 
processing with the calculation along a hyper-surface, which will lead to a min- 
imum. In fact, if the charge density is a unique functional for a given system, 
then the convergence criterium can be reformulated in terms of the difference 
of charge in two subsequent iteration cycles: if this difference is zero, then the 
charge density is the true groundstate charge density. Thus the problem of 
finding the solution of the non-linear eigenvalue problem and the groundstate 
charge density from: 



can be broken down in a number of discrete steps. Let us consider a general 
case; a system with a number M of ions, at the positions Ri of our coordinate 
system. The number of electrons N shall be: 



if the system is neutral and the ionic charge of atom i is given by Zi. If we 
neglect for the time being magnetic systems, then the potential v(r) is the sum 
of ionic and electronic contributions; it can be written as (we use atomic units 
in the rest of the paper) : 



Then the potential part of the eigenvalue equation is described by (we follow 
the convention in DFT and write the product of external potential and charge 



[-2/V 2 + v(r)p(r)] *(r) = etf (r) 



(38) 




(39) 




(40) 
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density as the effective potential w e //(r)): 

Veff (r) = v(r M r) = - f + J d *r>f^- (41) 

The eigenvalue problem then has the form: 

[-2/V 2 + v eff (r)] *(r) = ett(r) (42) 

A minimization process proceeds from one charge density to the next without 
the calculation of Kohn-Sham eigenstates. One may wonder that the parameter 
/, related to the kinetic energy density, still contains the volume. This is a 
consequence of the rigorous application of the energy principle, not only for 
the system as a whole, but for every single point. If the density of charge is 
increased, then the potential energy will also increase. This, in turn means 
that the kinetic energy density will be lower. Given a number of electrons in 
the system, all energy components depend on the volume. The volume must 
therefore also show up in the kinetic energy component. 

3.1 Initial charge density 

It is customary in DFT simulations to start a self-consistency cycle with an 
initial guess about the charge density distribution. In principle, this guess could 
involve any arbitrary choice, as long as the total number of electrons in the 
system is equal to N. However, most codes employ a superposition of atomic 
charges. As the previous section and the zero order approximation of the non- 
linear problem show, this is a suitable choice for the present problem. Since the 
charge density will ultimately be computed self-consistently, we may employ 
this approximation also for our initial charge distribution, setting po( r ) to: 

M Zi 

A>(r)=^5>cM(r-R,) (43) 

i— 1 a — 1 

where a indicates the atomic orbitals of a given atom at the position R^. From 
this initial charge density the potential v e f fp (r) is computed with: 

«.«.»W = -E0^ + /^^f^p (44) 

3.2 Self consistency cycle 

The initial potential is the input for the first iterative cycle to solve the ensuing 
linear density equation: 

[-2/V 2 + « e//;0 (r)]*(r) = e*(r) ( 45 ) 



15 



The problem is a standard problem in DFT, and a large number of algorithms 
exists for its solution. Note that it is only necessary to calculate the lowest 
eigenstate, defined by e%, and the corresponding real eigenvector pi(r). This 
is in marked contrast with today's methods, where an N electron system has 
to be solved for N eigenvalues and eigenvectors. This means, that the size of 
the problem depends only on the volume of the system. If we consider that 
the volume usually scales linearly with the number of atoms, then the ensuing 
numerical problem is also of the same size. Apart from the simplification of not 
having to solve the Schrodinger equation for single Kohn-Sham states the new 
method scales also much better than usual codes, which scale quadratically or 
cubic with N. It is a true order-TV method, or provides the optimum scaling 
with the size of the system. Today, order- N methods are based on a clever use of 
density matrices (see e.g |U), or they use atomic orbitals ^5], which makes them 
unsuitable for metallic systems. The suggested method has no limitation in this 
respect. Once the many-electron wavefunction is known, the charge density is 
given by: 

Pl (r) = **(r)* (r) (46) 

The density pi(r) is used to construct a new potential, u e //,i(r), by mixing 
the input density Po( r ) an d the output density p\(r) with the help of a suitable 
mixing algorithm, so that: 

Pi(t)=M\po(t),Pi(t)] (47) 

Again, this is standard procedure in today's DFT codes and a large number of 
mixing schemes (see e.g. [501121 EH) exists to guarantee numerical stability in 
the self-consistency cycle. The new potential v e ff,i(r) is then: 



f-f Ir - n t 



r — r 



And the second iteration cycle consists again of solving a linear eigenvalue prob- 
lem described by: 

[-2/V a (r)]*(r) = e*(r) (49) 

A solution of the problem leads to energy eigenvalue €2 and eigenvector ^2 (r) , 
which are used as an input for the new effective potential v e f f^(j). The iteration 
is repeated until convergency, usually defined either by the change of eigenvalues 
or the change of charge density, is achieved. 

3.3 Calculating the value of physical variables 

Once the groundstate charge density p g (r) has been computed the many-electron 
wavefunction can be determined by solving the eigenvalue problem: 

[-2/V 2 + v g (r) Pg (r)] *(r) = e*(r) (50) 
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The many-electron wavefunction may provide a direct way in the future to com- 
pute all physical quantities of interest. Since the many-electron wavefunction is 
the solution of the non-linear Hamiltonian, the density matrix: 

n(r,r') = #(r)#*(r') (51) 

can be used to calculate any physical quantity A of interest via: 

(A) = Tr [A* (r)#*(r')] = J d 3 rAp{r) (52) 

However, the physical quantities like kinetic energy and potential energy have 
been modified to obtain the general formulation. Similarly, other quantities in 
the Schrodinger equation, e.g. momenta, or electric field operators, will have to 
be adapted to the non-linear formulation. 

An alternative solution to the problem which allows to make use of the vast 
body of theory developed, will be to use only the many-electron groundstate 
charge density p g (r), and to construct the effective potential in the standard 
way by: 

w e //,#t(r) = v[p g (r)] (53) 

This effective potential then allows to construct single particle solutions in the 
standard manner: 

-^V 2 +«[ (09 (r)]|^(r)=6^(r) (54) 

These solutions can be used within the standard framework, considering that 
the groundstate charge density completely determines the physical properties 
of the system. The forces on ions in a system can be calculated in the usual 
manner. Since the Hamiltonian of the non-linear formulation is still hermitian, 
the Hellman-Feynman theorem still applies |23| : 



From this relation the forces on atomic nuclei can be calculated in the standard 
manner. These forces can then be used to optimize the geometry of the system. 



3.4 Summary 

We have shown in this section that the non-linear formulation of density func- 
tional theory allows to perform standard iteration procedures as used for at 
least three decades. It was pointed out that the main advantages of the new 
approach are: 

1. The iterations lead to the true many-electron charge density. 

2. The complexity of the calculation scales linearly with the number of atoms 
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Metal 


Type 


a [a.u.] 


V [a.u. 3 ] 


e/cell 


r s 


$[eV] 


Cs 


bee 


11.6 


1560.9 


2 


5.74 


1.81 


Li 


bee 


6.63 


291.4 


2 


3.26 


2.38 


Cu 


fee 


6.83 


318.8 


44 


1.20 


4.40 


W 


bee 


5.98 


214.0 


12 


1.62 


4.50 



Table 1: Lattice parameters, Wigner-Seitz radii of electrons, and workfunction 
of typical metals. The experimental workfunctions were taken from |24j . the 
lattice parameters from [2S1 



4 Numerical tests 

Lang and Kohn showed in the early Seventies 26 ;22j that the electron charge at 
the surface of a metal oscillates as it approaches the surface dipole. While stan- 
dard DFT is able to capture this behavior, the Friedel-oscillations, the Thomas- 
Fermi method fails quite dramatically. It seems therefore an ideal test case for 
non- linear DFT. The typical density of metals can be inferred from the lattice 
type, the lattice constant, and the number of valence electrons, assuming that 
the core covers only a small region of the unit cell. The Wigner-Seitz radii of 
typical metals, as well as their experimental workfunctions, are given in Table 
1. 

4.1 Setup 

We performed self-consistent simulations of the groundstate charge density for 
a metallic system, varying the density of the positive background in the jellium 
calculation from r s = 3 to r s = 1.5a. u. The metal was modeled by a slab of 
30(r s =1.5-2.0a.u.) to 40a.u.(r s = 2.5-3.0a.u.) thickness, corresponding to more 
than 15 layers of the metal. In standard DFT simulations this is sufficient to 
guarantee bulk properties in the center of the slab. In all calculations the film 
was embedded in two vacua, the vacuum range was 15-20 Bohr radii. This 
seems sufficient to mimic the surface potential barrier, it should also allow to 
estimate, whether the density decays exponentially into the vacuum range, as 
found in standard DFT simulations 26 , 3Jj , or inversely proportional with the 
distance from the jellium edge, as required by the asymptotic decay due to image 
potentials. The non-linear equation of the problem has the following form: 

<f(z) = E$(z) (56) 

Here, the eigenvalue E is the energy value of the many-electron electron charge. 
With a suitably chosen non-linear potential it will be equal to the eigenvalue 
of the many-electron Schrodinger equation. The dipole potential of the surface 
can be calculated in the standard way by integrating the contributions from a 



2dz 2 



V n l(z) 
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point deep inside the jellium into the vacuum region [251l27| : 

v el (z) = -4tt / dz 1 / dz" [p(z") - P+ (z")} , (57) 

J z J z' 

and by using the boundary condition v e i(z — > oo) = 0. The nonlinear potential 
was obtained by multiplying the dipole potential with V ws p{z), where V ws is the 
Wigner-Seitz volume: 

47T 

Vnl(z) = V e l{z)V ws p(z) V ws = —T s (58) 

The initial surface dipole was mimicked by a Fermi distibution function at both 
jellium edges. The initial density in the vacuum at a distance of 20a. u. was 
below 10 -17 in units of the jellium density. The initial nonlinear potential is 
nearly equal to the electrostatic potential inside the slab, the only significant 
differences occur at the jellium edges (see Fig. QJ. Conveniently, the nonlin- 
ear potential also has the same dimensions as the electrostatic potential. It 
starts deviating only, once the electron charge density in the initial setup differs 
substantially from the positive background. 



4.2 Iterations 

The nonlinear equation has been solved by numerical integration. We employed 
the Numerov integration scheme The value of ^ n +\ on an equally spaced 
z-grid is the convolution of the preceding ^-values and the second derivative of 
the equations, according to: 

_ 2*„ - *„_x + g (U n+1 + 10F n + 

1 h 'v n+1 ^ 59 J 

1 ' 12 

The precision of the result in one step is of the order h 5 , which renders the 
method quite accurate for one dimensional applications, h is the stepsize in the 
integration, the derivatives F n are given by: 

F„ = 2 {v nUn - E) (60) 

In case of the non-linear equation U n will be zero at every gridpoint, while V n is 
equal to 2 (v n i, n — E). The integration is performed twice: once starting from 
the left, integrating to the right, and once starting from the right and integrating 
left. The two functions are matched at z — 0. The parameter E, the eigenvalue 
of the system, determines whether they can actually be matched. The standard 
procedure in this case is to compute the derivatives of the two functions, ^i e ft 
and bright and to change the eigenvalue E until the two functions match. Nu- 
merical convergency with a simple linear mixing scheme required a very high 
number of gridpoints (8000 gridpoints over a length of 60-80a.u.) and very low 
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mixing parameters (0.001 to 0.002) in the simulation. The difference between 
input and output charge density is given by the integral: 



The charge density in the iterations was converged to a value of 0.01, which 
seemed sufficiently precise. The eigenvalue at this level of convergency is precise 
to about lmeV. The difference between input and output charge density in the 
iterations decreases exponentially, as shown in Fig. [21 

4.3 Density oscillations 

In the simulations we find, in every case, an oscillation of the charge density 
at the jellium edge. The wavelength of the oscillation depends on the jellium 
density. Contrary to Lang and Kohn 26, 27 we do not find that the oscillations 
arc damped for values of r s < 2.0a.u. The results of the simulations are shown 
in Fig. |2| We find the same increase of electron density at the jellium edge for 
every density value. The surplus charge is about 8% compared to the bulk value. 
The result suggests that within the non-linear framework a gradual removal of 
the peak and an approach of the Thomas- Fermi result in the limit of r s — > 2|26j 
does not occur. The standard DFT result for a uniform positive background 
cannot be extended beyond r s =3, as Lang and Kohn showed in their calculations 
[2111221 • This feature of the jellium model corresponds to a failure to account 
for the experimental surface energies in this range. At present, we have not yet 
developed a method to calculate surface energies from the many-electron charge 
densities alone. However, it seems quite possible that results remain reasonable 
also in this range, considering that the surface dipole seems independent of the 
background charge. Analyzing the details of the density oscillations into the 
bulk we find a shortening of the wavelength, as the density increases. This is 
mainly due to changes in the non-linear potential: as the potential becomes 
more negative, the eigenvalues do not keep up. The positive energy difference 
between the energy of the electron charge and the potential increases in this 
case, which leads to shorter wavelengths (see next sections). 

4.4 Potentials and eigenvalues 

The peak of the charge density at the jellium edge corresponds to a trough in the 
attractive potential. Compared to the initial potentials (see Fig. EJ the main 
changes are the deepening of the dip at the barrier and the potential oscillations 
decaying into the bulk (see Fig. |2J. Analyzing the actual process it has to be 
remembered that the many-electron wavefunction will be in phase throughout 
the film. As the two partial waves are matched at the center of the film, the 
actual wavelength is not completely arbitrary but reflects also the boundary 
conditions. From a physical point of view this is a somewhat unrealistic setup, 
because the waves will be reflected at the core of bulk atoms. The ensuing 




(61) 
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r s [a.u.] 


Eigenvalue [htr] 


$ [eV] 


Metal 


<S>exp[cV] 


3.0 


-0.055 


1.51 


Li 


2.38 


2.5 


-0.073 


1.98 


Mg 


3.64 


2.0 


-0.102 


2.80 


Al 


4.25 


1.5 


-0.160 


4.35 


W 


4.50 



Table 2: Eigenvalues and workfunctions for different densities. The workfunc- 
tions are in the range of 1.5 to 4.35eV, corresponding approximately to the 
range observed in experiments. Experimental values were taken from [23]. The 
deviation in the low density regime could be due to the neglect of the core region 
(see the text). 



multiple scattering effects should change the distribution of maxima and minima 
in a way not reproducible in a one-dimensional model. However, it is quite clear 
from the results that the many-electron wavefunction will have to be in-phase 
over a lengthscale of a few nanometer. The coherence can only be broken by 
elcctron-phonon interactions in a thermally activated environment. We would 
therefore conclude that long range effects are a necessary consequence of the 
present model. 

The eigenvalues obtained in the simulations are given in Table 2. They 
agree very well with experimental data, considering that the model is relatively 
simple. The only deviations are observed in the low density range r s > 3.0. 
This could be due either to neglecting the core volume in the calculation of 
the Wigner-Seitz radius, or to the fact that electron charge in a low density 
environment will experience the point- like and regular arrangement of the ionic 
cores much stronger than in a high density environment. In the first case, the 
actual Wigner-Seitz radius is smaller than the assumed value, in the second case 
the jellium model is no longer suitable. In experiments the measured values are 
in the range of 2.0 to 5.0eV. This corresponds to a WS radius in the range of 2.5 
to 1.3a.u. Comparing our results for the workfunction with full potential DFT 
simulations |29|. we note that the agreement with experiments for standard 
transition metals is improved. While in DFT simulations of bd metals the 
typical workfunctions are well above 5eV, in case of tungsten even above 6eV, 
they are still in the range of 4-5eV in the present model. We may conclude that 
the simulations reproduce the experimental values reasonably well, and that 
the problems occurring in jellium simulations of high density could be due to 
standard density functional theory itself. Physically speaking, one would expect 
that a high density of electron charge makes a continuous model more applicable 
than a low density environment. In standard DFT jellium simulations, one 
observes exactly the opposite trend [23 EZj ■ 



21 



4.5 Vacuum decay 

Another reason for obtaining lower workfunctions than in standard DFT could 
be the vacuum decay of electron density. The value for r s =3.0 at a distance of 
20a. u. from the surface in the fully converged system is about 10 _7 e/a.u. 3 . This 
is ten orders of magnitude higher than at the beginning of the self consistency 
loop. It is also three to four orders or magnitude higher than in standard DFT 
simulations. A logarithmic plot of the vacuum density reveals quite clearly, that 
the density does not decay exponentially in the long distance limit. Logarithmic 
plots of all systems are shown in Fig. EJ The curves deviate from an exponen- 
tial characteristic. The inset shows a fit of the density to two graphs: (i) An 
exponential with the decay constant equal to v 7 ^, and a function l/z. Both 
curves have the same value at the ultimate limit of the calculation, 15a. u. from 
the jellium edge. However, the exponential graph intersects the density curve 
at a discrete angle, while the l/z curve seems to match the asymptotic density 
decay. The match of the density to a l/z curve is limited to a small range 
of less than 5a. u. Even though the vacuum density is at least three orders of 
magnitude higher than in standard DFT simulations, the result is therefore not 
fully conclusive. We shall return to this point in future publications. 

4.6 Kinetic energy functionals 

The potentials, energy eigenvalues, and fc-values for different densities and two 
positions, at the center of the film, and at the peak of the surface dipole, are 
given in Table 3. The density, electrostatic potential, and nonlinear potential for 
r s = 1.5a.u. are shown in Fig. The first conclusion drawn from the numerical 
values is that the electrostatic potential is equal to the nonlinear potential at the 
center of the slab. This is understandable, if one considers that the density at 
this point is the initial density; the product of Wigner-Seitz radius and density is 
in this case equal to unity. Importantly, this means that there is no correction to 
the electrostatic potential at this point. The situation changes at the minimum 
of the surface dipole potential (z = zq), about 1.5a.u. from the jellium edge. 
Here, the charge density deviates from the initial value, and in this case the 
nonlinear potential will be affected. The difference, however is small and less 
than 10%, as the density peak does not exceed this surplus value. Concerning 
the eigenvalues, we find in all cases that they are equal to the potential at the 
center of the film. This means, that the kinetic energy density at the center 
of the film must be zero. A feature, which is only possible if the curvature of 
the wavefunctions at this point vanishes. The density oscillations at the jellium 
edge therefore decay into the bulk. Since the potential at the surface dipole is 
lower than the eigenvalue, the wavefunctions and also the densities will oscillate. 
Basing the oscillations, and the relation between fc-values and energies of a free 
electron gas in a potential v n i(zo), one would expect the A;- values to comply 
with the dispersion relation: 

h 2 = E~v nl {z ) (62) 
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r s 


Vel(0) 






Vnl(z ) 


E 


E - V n l(zo) 


k 


1.5 


-0.160 


-0.290 


-0.160 


-0.265 


-0.160 


0.105 


1.080 


2.0 


-0.102 


-0.190 


-0.102 


-0.174 


-0.102 


0.072 


0.910 


2.5 


-0.073 


-0.136 


-0.073 


-0.124 


-0.727 


0.051 


0.760 


3.0 


-0.056 


-0.104 


-0.056 


-0.095 


-0.556 


0.039 


0.660 



Table 3: Potentials at the center of the film z = 0, at the peak of the surface 
dipole z — Zo, eigenvalue E, difference between eigenvalue and attractive po- 
tential, and wavevector k. Note that the \k 2 is about five times larger than 
the difference between eigenvalue and attractive potential at z - The difference 
bears on the confinement of electrons, included in the kinetic energy functional 
(for details, see the text). 



As an inspection of the numerical data shows, this is not correct. In fact, the 
term on the left is about five times higher than the term on the right. This 
means, clearly, that the kinetic energy of the nonlinear equation is not the usual 
single-particle or Thomas-Fermi kinetic energy. The difference becomes clear if 
one remembers that the fundamental change in the Schrodingcr equation was 
dividing the kinetic single-electron term by the volume of the system. In a one- 
dimensional model, this corresponds to the lengthscale of the system. Here, we 
have to consider that the surface dipole is approximately 4-5a.u. wide: dividing 
the term on the left by 5 leads to the correct relation. Physically speaking, 
the change of the fundamental equation transformed the kinetic single-electron 
term into a kinetic energy functional. The decisive parameter of this functional 
is the confinement of the electrons in the field of the surface dipole. Regarding 
the energy eigenvalues it is quite clear that the energy density at the center of 
the slab is the same as the potential density. Approaching the jellium edge, the 
potential changes and the kinetic energy component changes accordingly: The 
system is indeed characterized by a constant energy density throughout. It is, 
therefore, a system in energetic equilibrium. The potential minimum between 
-7a. u and -3a. u. leads to an exponential decay of the wavefunction into the 
potential barrier and effectively screens the surface dipole from the bulk charge. 
This feature contains already the essence of a surface state, observed on many 
metal surfaces. 

4.7 Summary 

The method was tested for a metallic surface within a continuous background 
model, varying the density over a range from r s = 1.5 to r s = 3.0a.u. In this case 
we found, similarly to the results in standard DFT, that the density oscillates 
at the jellium edge, and that the wavelength of the oscillations depends on 
the jellium density. Wc compared the calculated eigenvalues with experimental 
values for the workfunctions and found good agreement, in particular in the 
high density range. The kinetic energy functional in the simulations includes 
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the confinement of electron charge. This has the effect that the energy density 
is constant throughout the system. Finally, we determined the decay of the 
density in the vacuum range. Here we find that the asymptotic behavior of 
the vacuum charge seems to be better described by a 1/z characteristic than 
an exponential function. In all simulations the only input parameters were the 
background density and the thickness of the slab. 
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Figure 1: (Color online) Setup for the calculation of metal surfaces. The jellium 
film is embedded in two vacuum ranges (top) . The initial electrostatic potential 
(full graphs) inside the slab is equal to the non-linear potential (dashed graphs) 
as defined in Eq. EH1 The potential varies with the density of the background 
(bottom). 



26 



t 1 1 1 1 r 




Iteration 

Figure 2: Convergency of the iterations. The difference between input and 
output charge density decreases close to exponentially to its final value below 
0.01. In the final iteration input and output charge density are identical (inset). 



27 




-15 -10 -5 5 
Distance from jellium edge [a.u.] 

Figure 3: (Color online) Density oscillations at the jellium edge. The density 
shows a characteristic peak at the boundary. The maximum value of the density 
is nearly equivalent in all simulations and about 8% higher than the density of 
positive charge. The oscillations decay into the bulk, the characteristic wave- 
length decreases with increasing density (see detail). 
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Figure 4: (Color online) Non-linear potential of the fully converged systems. 
The surface dipole corresponds to a minimum of the potential near the jellium 
edge. The characteristic wavelength A decreases with increasing density. 
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Figure 5: (Color online) Density decay in the vacuum region. The logarithmic 
plot shows a substantial deviation from an exponential behavior in the high 
distance regime. 



30 




-10 -5 5 

Distance from jellium edge [a.u.] 



Figure 6: (Color online) Density, electrostatic potential, and nonlinear potential 
for r s — 1.5a.u. The potential barrier around -4a. u. screens surface electrons 
from bulk electrons. 
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